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We study both static and transport properties of model quantum dots, employing density func- 
tional theory as well as (numerically) exact methods. For the lattice model under consideration 
the accuracy of the local-density approximation generally is poor. For weak interaction, however, 
accurate results are achieved within the optimized effective potential method, while for intermediate 
interaction strengths a method combining the exact diagonalization of small clusters with density 
functional theory is very successful. Results obtained from the latter approach yield very good 
agreement with density matrix renormalization group studies, where the full Hamiltonian consisting 
of the dot and the attached leads has to be diagonalized. Furthermore we address the question 
whether static density functional theory is able to predict the exact linear conductance through the 
dot correctly — with, in general, negative answer. 



I. INTRODUCTION 

Density functional theory (DFT) is an efficient tool for 
determining the electronic structure of solids. While orig- 
inally deve lop ed for continuum systems with Coulomb 
interaction}-^ DFT has also been applied to lattice mod- 
els, such as the Hubbard model or models of spinless 
fermionsP^ These lattice models often allow for exact 
solutions — either analytically or based on numerics — 
which hence can serve as benchmarks for assessing the 
quality of approximations. 

Very popular in solid-state applications is the local- 
density approximation (LDA) where the exchange- 
correlation energy of the inhomogeneous system un- 
der consideration is constructed via a local approxima- 
tion from the homogeneous electron systemP Recently 
a lattice version of LDA has been suggested for one- 
dimensional systems, where the underlying homogeneous 
system can be solved using the Bethe ansatz. For exam- 
ple, it has been demonstrated that the Bethe ansatz LDA 
describes well the low-frequency, long-wavelength excita- 
tions of the inte ract ing one-dimensional system, i.e., of a 
Luttinger liquidP^ 

On the other hand LDA often fails in correlated sys- 
tems and systematic improvements beyond the LDA are 
difficult. In this article we focus on a model of spinless 
fermions describing interacting electrons on a quantum 
dot. In a first step we compare the equilibrium proper- 
ties of the system, i.e., the number of particles on the dot 
as a function of the gate voltage obtained within differ- 
ent approximations for the exchange-correlation energy: 
the LDA and the optimized effective potential (OEP) ap- 
proach. Furthermore we suggest a novel method, where 
the exchange-correlation energy is obtained via the exact 
diagonalization of a small cluster that is composed of the 
strongly interacting region and a few additional sites. 

In the second step we compute the linear conductance 
through the dot. A general motivation for this study 
is recent progress in the field of molecular electronics, 
where DFT-based calculations are a standard tool to cal- 



culate electrical conductancesJ^HUl however the concep- 
tional limitations of the approach are not very well un- 
derstood yet. More specifically we were motivated by 
the model studies in Refs. [13] and [14]. Mera et al™ 
stressed that static DFT reproduces the conductance of 
an interacting system correctly if there exists a Friedel 
sum rule that relates the conductance with the equilib- 
rium density. Schmitteckert and Evera^ compared con- 
ductances obtained from a density matrix renormaliza- 
tion group calculation with those obtained within static 
DFT. Close to resonances both conductances were in very 
good agreement, while off-resonance there was a consider- 
able discrepancy. This discrepancy is due to a exchange- 
correlation contribution to the voltage difference between 
the two reservoirspl U xc . One of the questions we will 
address is whether U xc depends on the distance between 
the interacting region and the reservoirs. 

In the following section we will introduce the model 
under investigation. In Sec. Ill, devoted to static density 
functional theory, we discuss the approximations used 
to obtain the exchange-correlation energy. Section IV 
is devoted to transport: We rederive the Meir-Wingreen 
formula for the conductance starting from the dynami- 
cal density-density response function, and we apply the 
formula to calculate the conductance for our model. The 
final section contains a summary as well as our conclu- 
sions. 



II. THE MODEL 



We study a model where a chain of TV lattice sites is 
coupled to two reservoirs 

H = H L +H LC + Hc+ H CR + Hr. (I ) 
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FIG. 1: Schematic view of the model under consideration: A 
quantum dot (full circles) attached to left and right leads, each 
of them consisting of m non-interacting sites (open circles) 
and a reservoir described by a continuum of states (shaded 
regions). The electrons on the dot (site m + 1 to m + 5) 
interact and their potential energy can be tuned by a gate 
voltage. 



The Hamiltonian of the chain reads 
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where c/~ and di are fermion creation and annihilation 
operators, hi = c^q counts the fermions on lattice site /. 
The hopping matrix elements are chosen such that 

the system resembles a quantum dot that is weakly cou- 
pled to left and right leads, cf. Fig. [T] and are explicitly 
given by 




. , m — 1, m + 6, . 
and I = m + 5 
1, . . . , m + 4 



,iV-l 
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In the following we choose t' = 0.2t and t<jot = 0.5t. 
The interaction strength and the potentials are constant 
within the quantum dot, V^+i = V and vi — u ga tc, and 
zero outside. The reservoirs, chosen to be non-interacting 
fermions, are described by 
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E 6k ^kL/R^kL/R^ 



(4) 



and finally the coupling between the chain and the left 
reservoir reads 
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LC 



C kL C l + C l C kL 



(5) 



and analogously for Hqr- The hopping parameters t k 
are fine-tuned such that an electron at the Fermi energy 
is not back-scattered at the interface between the chain 
and the reservoirs = ^t/irN(ek) where N(ek) is the 
density of states in the reservoirs. Up to this point our 
model is the same as the model studied by Schmitteckert 
and Evers [T3]- In contrast to [T3], however, we replace 
the discrete levels in the reservoirs by a continuum of a 
wide and flat band, which then leads to 

T L/R = 27rY,\tk\ 2 S(e F ~e k ) (6) 

k 

= 2ir\t k \ 2 N(e k )\ =2t. (7) 



Notice that in the continuum limit the density of states 
in the reservoirs goes to infinity, and thus — in order to 
keep T L / R constant — the coupling strength between the 
states in the reservoirs and the chain goes to zero. 



III. STATIC DENSITY FUNCTIONAL THEORY 

The lattice version of DFT relies on the fact that there 
is a one-to-one correspondence between local potentials 
{v^ and the ground-state expectation values of the site 
occupations {n,}. Therefore it is — in principle — possible 
to express all quantities that can be obtained from the 
ground-state wave function as a function of the densities. 
The site occupations as a function of the potentials can 
be found from derivatives of the ground-state energy with 
respect to the local potential 
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In order to determine the potentials from the densities it 
is convenient to define the function 



F ({ n i}) = mm*_ 



;*|f + U|*> 



(9) 



where \& — ¥ {n^} indicates that the minimization is con- 
strained to such wave functions W that yield the given site 
occupations {rii}. Here T and V are the kinetic and the 
interaction parts of the Hamiltonian, respectively. The 
ground-state energy is obtained by minimizing the func- 
tion 



E({n i }) = F({n i }) + J2 



v l n i 



(10) 



with respect to rij. When we minimize E under the con- 
straint of a constant particle number we obtain the poten- 
tial up to an additive constant (the Lagrange multiplier) : 



dF 
dm 



+ A. 



(11) 



A major step towards the practical implementation of 
DFT is to employ a non-interacting auxiliary Hamilto- 
nian H s (Kohn-Sham Hamiltonian) in order to calculate 
the density profile, 



(12) 



where the potentials i?| have to be chosen such that in the 
ground-state of H s the site occupations are the same 
as in the interacting model. In analogy to the interact- 
ing system, the ground-state energy of the Kohn-Sham 
system is found by minimizing 



E*({n i }) = F s ({n i }) + Yv!n i 



(13) 
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Combining ( 10 1 and ( 13 ) yields 



with the Fock-like exchange energy 



E({ ni }) =E s ({n l })+E }ixc an l })+J2(v t -v!)n l , (14) 



(19) 



with the Hartree-exchange-correlation energy defined by We ^ calculate the ground state expectation values 



E 



HXC 



({n i }) = F({n i })-F*({n i }). 



(15) 



The condition that both E and E s are minimal for the 
same set of site occupations ni requires that 



dm 



(16) 



Up to this point no approximations have been em- 
ployed. However to determine £ HXC at a given density 
exactly is as demanding as finding the ground state en- 
ergy for a given potential. The hope is that there exist 
good approximations for £? HXC that are accessible with 
low numerical cost but still allow good estimates for the 
ground-state energy and density. Here and in the fol- 
lowing we will compare three different approximations: 
the local-density approximation, the optimized effective 
potential (so-called exact exchange) approximation, and 
finally a method based on the exact diagonalization of 
small clusters. 



A. Local density approximation 

In the LDA one writes E KXC as the sum of the (non- 
local) Hartree energy plus an exchange-correlation energy 
which depends only on the local density, 
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5> x cK). (17) 



The local exchange-correlation energy exc('T-) is deter- 
mined from the ground-state energy density of a ho- 
mogeneous system at the same density. For the one- 
dimensional lattice models this quantity can be calcu- 
lated using the Bethe ansatz, see Refs. |6I8| . Notice that 
in the Hamiltonian we study the interaction strength de- 
pends on position, and an ambiguity arises how to deter- 
mine the exchange-correlation potential for those sites 
which interact only with one neighbor. For simplicity we 
used in our numerical implementation of the LDA the 
same function e xc (n) for all interacting lattice sites, i.e., 
number m + 1 to m + 5. 



B. Optimized effective potential 

In the OEP approach the Hartree-exchange-correlation 
energy is 

£qep ({»*}) = Vj2n t n l+1 +E x ({ ni }) (18) 



Cji i) etc. for a non-interacting system coupled to 
reservoirs using the Green's function technique, see sec- 
tion IV. The Kohn-Sham equations ( 16 ) are most conve- 



niently solved by iteration. Starting with an initial guess 
for the potentials vf we calculate the corresponding site 
occupations and the Hartree-exchange-correlation en- 
ergy E HXC . Since in the OEP approach £; HXC depends 
only implicitly on m we rewrite Eq. ( 16 ) using 



d£ HXC 
dm 



E 



dE nxc dv s 
dv? drij 



(20) 



where the derivatives of £' HXC with respect to the v? are 
calculated numerically and dvj /drii is obtained by ma- 
trix inversion from drii/dv^. Finally we obtain a new 
set of Kohn-Sham potentials vf. The whole procedure is 
repeated until convergence, i.e., until the difference be- 
tween old and new potentials is smaller than some given 
cutoff. 

The fact that we use a continuum of states to de- 



scribe the reservoirs simplifies the task to solve Eq. ( 16 ) 



considerably: Since the states in the reservoirs are only 
infinitcsimally weakly coupled to the dot, the Hartree- 
exchange-correlation potential in the reservoirs disap- 
pears, so that the number of potentials vf to be deter- 
mined self-consistently equals the chain length N. 



C. Exact diagonalization 

In our model Hamiltonian electrons interact only in a 
spatially confined region. In the non-interacting regions 
we find numerically (e.g., within the OEP approach, see 
Fig. [6] below; compare also Ref. [13] for density ma- 
trix renormalization group (DMRG) results) only small 
exchange-correlation potentials. This finding motivated 
us to use the exchange-correlation energy of a small clus- 
ter consisting of the interacting region plus a small num- 
ber of non-interacting sites as an approximation for the 
exchange-correlation energy of the system attached to 
reservoirs: 



E 



HXC 
ED 



({ni}) = F^ttm}) - F£ D ({ ni }) (21) 



where F ED ({rii}) and F ED ({rii}) are exact on the small 
cluster and can be obtained by numerical diagonalization. 

In this approach we have to fine-tune the local po- 
tentials of three different Hamiltonians such that all the 
three yield the same local densities: (i) a cluster of in- 
teracting electrons with potentials Ui, (ii) a cluster of 
non- interacting electrons with potentials uf, and (iii) 
the Kohn-Sham Hamiltonian of the extended quantum 
dot attached to reservoirs with vf — Vi + i> HXC where 
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FIG. 2: Particle number on the dot, iVdot, as function of the 
gate voltage, u ga t c , for V/t = 0.25 obtained within density 
functional theory with three different approximations for the 
exchange-correlation energy: the local-density approximation 
(LDA), optimized effective potential (OEP) and a method 
based on the exact diagonalization of short chains (ED). Here 
the chains have a length of nine lattice sites. 
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FIG. 4: iVdot as function of u ga to for V/t = 0.5. Solving the 
OEP equations iteratively as a function of gate voltage two 
solutions are found close to « ga tc = 0. 
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FIG. 5: iVdot as function of « ga te for V/t = 2. For comparison 
we also include the DMRG results of 1131. 



FIG. 3: iVdot as function of w ga t c around to the central step 
for V/t = 0.25; for this interaction strength OEP is still close 
to the exact result. 



vf xc = Ui — uf on the cluster sites and zero outside. 
Again, in a practical scheme the fine-tuning procedure is 
performed by iteration. The task to determine the poten- 
tials Ui that correspond to a given set of site occupations 
rii is nontrivial for an interacting system and limits the 
cluster size in our approach to approximately 12 to 14 
sites. 
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D. Results 

In the following we focus on the particle number in 
the interacting region, N^ot, as function of the gate volt- 
age, Vgato, comparing results for the three aforementioned 
approaches LDA, OEP, and ED, respectively. The data 
are obtained for a chain of nine sites, i.e., the five-site 



FIG. 6: Hartree-exchange-correlation potential up xc in OEP 
for a long chain with m = 20 non-interacting sites attached 
on each side of the five-site quantum dot with parameters 
Vgntc/t — 0.5 and V/t — 0.5. In the figure, up xc is shown for 
the non-interacting sites I — 1, ... ,m and for the whole chain 
in the inset. 
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quantum dot plus two non-interacting sites on each side 
of it. Fig. [2] shows iVdot for weak interaction strength, 
V/t = 0.25. The three curves nearly coincide with excep- 
tion of the regions close to the steps, in particular around 
V gatc = 0. Here, as pointed out in Fig. [31 the step ap- 
pears steeper in OEP and flatter in LDA compared to 
ED. 

This trend continues at stronger interaction, V/t = 0.5, 
where in OEP the particle number even jumps close to 
Vgatc = with a small hysteresis region of two stable 
solutions, while the LDA step flattens out even more, as 
displayed in Fig. |4| 

In the strong coupling regime (Fig. [5| where V/t = 
2, comparison with the exact densities obtained from 
DMRG calculations^ shows that LDA and OEP fail 
completely, while the ED results agree reasonably well 
with the exact data. Note however that the DMRG data 
are for a five-site quantum dot with two non-interacting 
sites attached on the right and three on the left, and 
that the reservoirs are described by a finite set of the or- 
der of 100 discrete levels instead of a continuum. From 
the above observations we conclude that in DFT calcula- 
tions for lattice models LDA and OEP results are reliable 
in the weak interaction regime. In particular, it can be 
shown that the OEP density profile is exact to linear or- 
der in the interaction strength V. On the other hand, for 
strongly correlated systems more sophisticated methods 
like the ED cluster approach are required, even for static 
properties. 

As already mentioned in Sec. III.C, far from the in- 
teracting region the potential i> ; HXC becomes very small. 
This is explicitly demonstrated in Fig. [6j In the leads 
the Hartree-exchange-correlation potential, w ; HXC , within 
OEP, is found to be about three orders of magnitude 
smaller than in the interacting region (see inset). 



IV. TRANSPORT 



particle number in the reservoir, 



I = e N L =eY^ 



rikL, 



(22) 



i.e., = —ieuNi,(uj), where — e is the electron charge. 
The frequency dependent variation in the particle num- 
ber Nl (w) appears as a response to a perturbation in the 
(single particle) Hamiltonian of the form 



(23) 



where Sv^ — <5v° xt + 5v^ xc is the sum of an external 
potential and the induced Hartree-exchange-correlation 
potential. The summation a includes both reservoir, 
a = kL, kR and chain degrees of freedom, a = I. The 
variation of the density at site j3 is then 



Sri/3 (u) 
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(0,a)Sv s a (u;)g e (a,/3), (24) 



where Q € (a,0) is the (zero temperature) Green's func- 
tion of the single-particle Hamiltonian. It is useful to 
distinguish the Green's function of the reservoirs from 
the Green's function of the chain, and in the following 
we will use the symbols g e (kL) for the (left) reservoir 
Green's function and G e (l,l') for a chain Green's func- 
tion. The latter is given by 

G e (l, I') = G° e {l, + G °^' m)^ e {m)G e (m, I'), (25) 

m 

were G®(l,l') is the bare Green's function, i.e., the one 
for tk — 0, and the self-energy S e (m) appears due to the 
coupling to the reservoirs. For our model Hamiltonian 
S e (m) is non-zero only on the first and last site of the 
chain, m — 1, N. The explicit expression for the first site 
is 



x e (i) = J2\tk\ 2 g°AkL), 



(26) 



DFT as presented in the previous sections is a ground 
state theory. However generalizations are available which 
allow to calculate densities at finite temperature and 
under non-equilibrium conditions^ and thus — via the 
continuity equation — the current through the quantum 
dot. The goal of this section is to calculate the DC- 
conductance of the quantum dot when a small voltage dif- 
ference is applied. We will extract the DC-conductance 
from a calculation of the dynamic density response func- 
tion, and the connection to the standard Meir-Wingreen 
formula 16 will be made. 



A. Conductance from density response 

We start with the current flowing from the left reservoir 
into the dot, which is given by the time derivative of the 



where 
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(27) 



is the bare Green's function of state k in the left lead. 
For energies close to the chemical potential and for our 
special choice of the couplings t k the self-energy then 
assumes the value £ e (l) = — itsga(e — /z). The Green's 
function for the states in the left reservoir finally is 

9t (kL) = 9 ° e (kL) + g° e (kL) t k G e (l, 1) t k g° e (kL). (28) 

The variation of the particle number in the reservoirs can 
be represented graphically as 



5N L = 




b) 



(29) 



G 



where the full line denotes a reservoir Green's function, 
g®(kL) or g®(kR), the broken line is the chain Green's 
function G e (l, I') and the cross corresponds to a hopping 
process between reservoir and chain. In order to calculate 
the DC conductance of the system we have to evaluate 
these diagrams for small but finite frequency, and we have 
to identify contributions that diverge as 1/ui as w goes to 
zero. Such divergences are found in diagrams c), d) and 
c). In all three cases the diverging contribution arises 
from the region in the e- integration where e < fi < e + u>. 
For diagram c), for instance, the relevant contribution is 
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G e + U (l, 1) 



(30) 
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In the next step the fc-summation in the second line is re- 
placed by an integral. Clearly the dominant contribution 
to the /c-integration comes from a small region around 
the Fermi energy. Assuming that neither the potential 
Svf. L nor the coupling t k are singular around the Fermi 
momentum we find J2k(---) = ^l8vI l /lj 2 , where the 
potential has to be evaluated at the Fermi energy. The 
variation of the particle number in the left reservoir hence 
is 



«v2 = ^ 

L 2ttw 



r £ G*(i,i). 



(31) 



The retarded chain Green's function at the che mica l po- 
tential, G^(l, 1), appears since e + oj > fi in Eq. (30) and 
we consider the limit uj — > 0. 

Using similar arguments for the diagrams d) and e), 
we find 
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[r L G*(l,N)r R G$(N,l)]. (33) 
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Using finally the relation 

g*(i,i)-g£(i,i) = -iG*(i,i)r L G£(i,i 



-iGK(l,N)r R G*(N,l)(M) 



the complete singular contribution to the density in the 
reservoir reads 



SN L = (5v s kR -5v s kL )—r L GK(l,N)r R G$(N,l). (35) 

This enables us to write the current as the product of a 
conductance G and a voltage U tot , 



I = -ieojSN L = GU to \ 



(36) 



where the expression for the conductance agrees with 
the standard Meir-Wingreen formula for non-interacting 
electrons, 



G = —r L G*(l, L)T R G$(L, 1), (37) 



and the voltage is the sum of an externally applied volt- 
age and an exchange-correlation contribution, U tot = 
L7 cxt + f/ xc , with 



rrCXt r cxt e cxt 

eu — ovu ' 



kR -5v™, eU™ = 5v^-Sv^. (38) 



Notice that only the exchange-correlation potentials in 
the reservoirs but not in the chain contribute to the total 
voltage. Our result is thus consistent with [T7] where 
it has been shown that the exact current can be ex- 
pressed in terms of a Landauer-type formula^ in which 
the electro-chemical potential of the leads is shifted by 
the volt age- induced variation of the exchange correlation 
potential, and with the statements of Refs. |14ll9j that 
static DFT gives the exact linear-response conductance 
provided that the dynamic exchange-correlation poten- 
tial vanishes deep inside the leads. 

Notice that U xc is a purely dynamic effect that cannot 
be captured by any adiabatic approximation. However, 
L7 XC can be assessed via reverse engineering: If we know 
the exact static densities and the exact conductance, then 
the ratio of the exact conductance and the DFT conduc- 
tance, Eq. ( [37] ), is equal to the ratio of the total and the 
external voltage, G c ™ ct /G = U tot /U cxt . 



B. Results 

Figure [7] shows the conductance as a function of the 
gate voltage in the case of a weak interaction V — 0.25£. 
One observes five resonances at the gate voltages where 
the particle number on the dot changes, compare Fig. [2] 
LDA overestimates the width of the resonances just as 
it overestimates the width of the steps in the particle 
number. We also include OEP and Hartree-Fock (HF) 
results in the figure. For weak interaction both methods 
predict identical charge densities, for the present inter- 
action strength the difference in particle number on the 
dot is for the two methods less than 7 • 10~ 3 . Also the 
conductances are close to each other. 

In Fig. [H] we show the region close to zero gate voltage 
in more detail. Near the resonance HF and OEP are al- 
most indistinguishable, however far from the resonance a 
significantly different conductance is found. This differ- 
ence is due to the exchange-correlation contribution to 
the voltage, U xc . To substantiate this point we analyze 
the correction to the conductance to first order in the 
interaction strength V. In this case HF yields the exact 
conductance. 

However, as demonstrated in Fig. [9] the conductances 
obtained by OEP and HF differ — although the densities 
are identical. Figure [9] has been obtained for a chain 
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FIG. 7: Conductance G as function of v ga .te for V/t = 0.25. 
The OEP and LDA curves are obtained within static density 
functional theory as explained in the text. "HF" corresponds 
to a self-consistent Hartree-Fock calculation. 




0.01 



-0.3 -0.2 -0.1 0.1 0.2 0.3 

Wgate /t 

FIG. 8: Conductance G as function of u ga te close to the central 
peak for V/t = 0.25. 
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FIG. 9: Linear contribution G^' of the expansion of the con- 
ductance G = G (0) + G w V/t + G (2) (V/t) 2 + ... in powers 
of the interaction strength. The difference between the OEP 
result and the exact HF conductance, AG' 1 ' = Gqe P — G^ F , 
is plotted as function of u ga t c . 




FIG. 10: Conductance G as function of v ga te for V/t = 2; the 
DMRG data are taken from Ref. 1131. 



length N — 9, but we have checked that even for much 
longer chains (up to N = 25) there is no visible change in 
the results. This means that U xc remains nonzero even 
far from the interacting region. 

Figure [10] shows the conductance for the relatively 
strong interaction strength V = 2t, comparing LDA, ED 
and the numerically exact conductance obtained with the 
DMRG in Ref. [13 . Although the position of the reso- 
nances is not bad in LDA the method predicts conduc- 
tances that differ several orders of magnitude from the 
exact results. The ED results are considerably better. 
Near the resonances ED predicts conductances that are 
close to the DMRG values. 



V. SUMMARY AND CONCLUSIONS 

We studied the ground state density profile and the 
conductance of a model quantum dot comparing DFT 
and exact results. The electron density in the ground 



state can be obtained reliably using non-local exchange- 
correlation potentials. While for weak interaction the 
OEP approach gives good results, in the case of interme- 
diate or strong interaction strengths a non-local potential 
extracted from the exact diagonalization of small clusters 
(DFT+ED) works well. 

For the conductance our results are not so clear-cut. In 
our simple model we find five well separated resonances 
as a function of the gate voltage. Static DFT reproduces 
very well the position and the width of these resonances. 
However, in contrast to model cases where the Friedel 
sum rule allows to express the conductance in terms of 
the equilibrium charge density, compare Ref. |14| . in our 
model system such a sum rule is not valid. As a con- 
sequence the DFT conductance is not exact, and we 
find pronounced deviations in the conductance valleys 
between the resonances. 

As the origin of these discrepancies we have identi- 
fied a dynamical exchange-correlation correction to the 



applied voltage, which is non-zero even when the reser- 
voirs are far from the interacting region. We believe that 
this finding is related to "ultra-non-locality" which is an 
inherent problem in time-dependent density functional 



theorjt™. 
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